# Load Medland survey collection data
medland_survey2014_2017 <- read_csv("data/medland_survey2014_2017.csv", locale = locale(encoding = "850"))
Rows: 1501 Columns: 28── Column specification ───────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): study.area, sector, subsector, visibility, time.date
dbl (23): ID, zone, collection, xcoord, ycoord, area.sqm, undiag.lithics, burins, notch.dent, MP.to...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
medland_survey2014_2017 %>%
dplyr::filter(area.sqm>0 & total.lithics>0 & !is.na(visibility)) %>%
mutate(lithic.density=total.lithics/area.sqm) %>%
ggplot(aes(x=lithic.density), xlim=.02) +
geom_histogram(binwidth = .001) +
scale_y_log10() +
scale_x_continuous(limits = c(0,0.01), breaks = c(0,.002, .004, .006, .008)) +
labs(title = "Surface Visibility and Artifact Recovery",
x="lithic artifacts / km^2",
y='count of patches') +
facet_grid(factor(visibility)~study.area) +
theme_bw(base_size = 20) +
theme(axis.title.x = element_markdown())
cat("\nCanal de Navarrés survey area\n")
Canal de Navarrés survey area
with(medland_survey2014_2017 %>%
dplyr::filter(total.lithics>0 & !is.na(visibility) & study.area == "Canal de Navarrés") %>%
mutate(lithic.density=total.lithics/area.sqm),
aov(lithic.density~visibility)) %>% summary()
Df Sum Sq Mean Sq F value Pr(>F)
visibility 2 0.0000022 1.104e-06 0.451 0.637
Residuals 194 0.0004746 2.446e-06
cat("\nHoya de Buñol survey area\n")
Hoya de Buñol survey area
with(medland_survey2014_2017 %>%
dplyr::filter(total.lithics>0 & !is.na(visibility) & study.area == "Hoya de Buñol") %>%
mutate(lithic.density=total.lithics/area.sqm),
aov(lithic.density~visibility)) %>% summary()
Df Sum Sq Mean Sq F value Pr(>F)
visibility 2 1.220e-06 6.095e-07 0.4 0.671
Residuals 98 1.493e-04 1.523e-06
cat("\nCocina/Catadau survey area\n")
Cocina/Catadau survey area
with(medland_survey2014_2017 %>%
dplyr::filter(total.lithics>0 & !is.na(visibility) & study.area == "Cocina/Catadau") %>%
mutate(lithic.density=total.lithics/area.sqm),
aov(lithic.density~visibility)) %>% summary()
Df Sum Sq Mean Sq F value Pr(>F)
visibility 2 0.0001894 9.470e-05 2.662 0.0758 .
Residuals 83 0.0029531 3.558e-05
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Modify training data based on ML and Bayesian testing:
# Merge ENEOL and MNEOL
# Remove undiagnostic lithics for age estimates
# Sort factor levels chronologically
training_data_E_Iberia <- training_data_E_Iberia %>%
select(-undiag.lithics, -total.lithics, -citation) %>%
mutate(period = replace(period, period == "MNEOL", "ENEOL"),
period = factor(period,
levels = c("MP", "UP", "EPI", "MESO", "ENEOL", "LNEOL")))
# filter out assemblages with only undiagnostic lithics
medland_survey2014_2017_lithics <- medland_survey2014_2017 %>%
dplyr::filter(undiag.lithics < total.lithics) %>%
select(ID, c(13:27))
# make separate table of assemblage ID and provenience
medland_survey2014_2017_info <- medland_survey2014_2017 %>%
dplyr::filter(undiag.lithics < total.lithics) %>%
mutate(assemblage = paste(study.area, "-", zone, "-", sector, "-", subsector, sep = "")) %>%
select(ID, study.area, assemblage)
# calculate lithic density for each collection patch
medland_survey2014_2017_density <- medland_survey2014_2017 %>%
dplyr::filter(undiag.lithics < total.lithics) %>%
mutate(assemblage = paste(study.area, "-", zone, "-", sector, "-", subsector, sep = ""),
density.km2 = total.lithics*1000/area.sqm)%>%
select(ID, total.lithics, area.sqm,density.km2)
# Partition into training and hold out test / validation sample
set.seed(456) ## if we want to make it completely reproducible
vl.split <- training_data_E_Iberia %>%
rsample::initial_split(., prop=.75)
vl.train <- rsample::training(vl.split)
vl.test <- rsample::testing(vl.split)
# save ID data for later analysis
vl.test.id <- vl.test %>%
select(ID, period)
10 folds using all the training data
set.seed(456)
vf10.all <- vfold_cv(training_data_E_Iberia %>% select(-ID),v=10)
10 folds using the 75% training data split
set.seed(456)
vf10.train <- vfold_cv(vl.train %>% select(-ID),v=10)
valencia.lithics.rf.mod <-
rand_forest(trees=500) %>%
set_engine("ranger", importance = "impurity") %>%
set_mode("classification")
print(valencia.lithics.rf.mod)
Random Forest Model Specification (classification)
Main Arguments:
trees = 500
Engine-Specific Arguments:
importance = impurity
Computational engine: ranger
valencia.lithics.rf.fit <-
valencia.lithics.rf.mod %>%
fit(as.factor(period) ~ ., data = vl.train[,2:ncol(vl.train)])
print(valencia.lithics.rf.fit)
parsnip model object
Ranger result
Call:
ranger::ranger(x = maybe_data_frame(x), y = y, num.trees = ~500, importance = ~"impurity", num.threads = 1, verbose = FALSE, seed = sample.int(10^5, 1), probability = TRUE)
Type: Probability estimation
Number of trees: 500
Sample size: 96
Number of independent variables: 15
Mtry: 3
Target node size: 10
Variable importance mode: impurity
Splitrule: gini
OOB prediction error (Brier s.): 0.2416893
valencia.lithics.rf.fit %>%
extract_fit_engine() %>%
vip(aesthetics = list(color = "black", fill = "#26ACB5"), num_features = 15) + theme_minimal()
valencia.lithics.rf.predicted <-
valencia.lithics.rf.fit %>%
predict(vl.test) %>%
bind_cols(vl.test.id[1:2], ., predict(valencia.lithics.rf.fit, vl.test, type="prob")) %>%
rename(predicted.age = .pred_class,
MP=.pred_MP,
UP=.pred_UP,
EPI=.pred_EPI,
MESO=.pred_MESO,
ENEOL=.pred_ENEOL,
LNEOL=.pred_LNEOL,
true.age = period) %>%
mutate(true.age = factor(true.age,
levels=c("MP", "UP", "EPI", "MESO", "ENEOL", "LNEOL")),
predicted.age = factor(predicted.age,
levels=c("MP", "UP", "EPI", "MESO", "ENEOL", "LNEOL")))
print(valencia.lithics.rf.predicted)
valencia.lithics.rf.predicted %>%
pivot_longer(cols = 4:ncol(valencia.lithics.rf.predicted),
names_to = "period", values_to = "probability") %>%
mutate(period = factor(period,
levels = c("MP", "UP", "EPI", "MESO", "ENEOL", "LNEOL"))) %>%
ggplot(aes(x=period, y=probability)) +
geom_line(group=1) +
geom_vline(aes(xintercept = true.age), color="red", size=2, alpha=.5) +
geom_vline(aes(xintercept = predicted.age), color="blue", size=0.8) +
labs(title="Random Forest Predictions for Each Assemblage",
subtitle="black line indicates probability, blue line indicates predicted age, & red line indicates known age",
x="time period",
y="probability of predicted time period") +
facet_wrap(vars(ID), ncol = 7) +
theme_bw(base_size = 20) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
strip.text.x = element_text(size = 14))
library(caret)
with(valencia.lithics.rf.predicted,
caret::confusionMatrix(true.age, predicted.age))
Confusion Matrix and Statistics
Reference
Prediction MP UP EPI MESO ENEOL LNEOL
MP 13 0 0 0 0 0
UP 0 2 0 0 2 0
EPI 0 1 1 1 0 0
MESO 0 0 0 3 1 0
ENEOL 0 0 0 1 5 0
LNEOL 0 0 0 0 2 1
Overall Statistics
Accuracy : 0.7576
95% CI : (0.5774, 0.8891)
No Information Rate : 0.3939
P-Value [Acc > NIR] : 2.407e-05
Kappa : 0.6788
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: MP Class: UP Class: EPI Class: MESO Class: ENEOL Class: LNEOL
Sensitivity 1.0000 0.66667 1.00000 0.60000 0.5000 1.00000
Specificity 1.0000 0.93333 0.93750 0.96429 0.9565 0.93750
Pos Pred Value 1.0000 0.50000 0.33333 0.75000 0.8333 0.33333
Neg Pred Value 1.0000 0.96552 1.00000 0.93103 0.8148 1.00000
Prevalence 0.3939 0.09091 0.03030 0.15152 0.3030 0.03030
Detection Rate 0.3939 0.06061 0.03030 0.09091 0.1515 0.03030
Detection Prevalence 0.3939 0.12121 0.09091 0.12121 0.1818 0.09091
Balanced Accuracy 1.0000 0.80000 0.96875 0.78214 0.7283 0.96875
#Create workflow step
valencia.lithics.rf.wf <-
workflow() %>%
add_model(valencia.lithics.rf.mod) %>%
add_formula(period ~ .) #The predictor is contained in add_formula method
set.seed(456) # For reproducibility
valencia.lithics.rf.xv.fit <-
valencia.lithics.rf.wf %>%
fit_resamples(vf10.all,
control=control_resamples(save_pred = TRUE))
! Fold01: internal: No observations were detected in `truth` for level(s): 'UP', 'LNEOL'
Computation will...
! Fold02: internal: No observations were detected in `truth` for level(s): 'EPI', 'LNEOL'
Computation wil...
! Fold03: internal: No observations were detected in `truth` for level(s): 'EPI', 'LNEOL'
Computation wil...
! Fold04: internal: No observations were detected in `truth` for level(s): 'UP'
Computation will proceed ...
! Fold06: internal: No observations were detected in `truth` for level(s): 'EPI', 'MESO'
Computation will...
! Fold09: internal: No observations were detected in `truth` for level(s): 'LNEOL'
Computation will proce...
! Fold10: internal: No observations were detected in `truth` for level(s): 'EPI', 'LNEOL'
Computation wil...
print(valencia.lithics.rf.xv.fit)
# Resampling results
# 10-fold cross-validation
There were issues with some computations:
- Warning(s) x3: No observations were detected in `truth` for level(s): 'EPI', 'LNEOL' Computation wi... - Warning(s) x1: No observations were detected in `truth` for level(s): 'EPI', 'LNEOL' Computation wi... - Warning(s) x1: No observations were detected in `truth` for level(s): 'EPI', 'LNEOL' Computation wi... - Warning(s) x1: No observations were detected in `truth` for level(s): 'EPI', 'LNEOL' Computation wi... - Warning(s) x1: No observations were detected in `truth` for level(s): 'EPI', 'LNEOL' Computation wi...
Use `collect_notes(object)` for more information.
# Collect the metrics using another model with cross-validation
valencia.lithics.rf.xv.meanpreds <- tune::collect_metrics(valencia.lithics.rf.xv.fit)
print(valencia.lithics.rf.xv.meanpreds)
valencia.lithics.rf.mod %>%
fit(as.factor(period) ~ ., data = training_data_E_Iberia[,2:ncol(training_data_E_Iberia)]) %>%
extract_fit_engine() %>%
vip(aesthetics = list(color = "black", fill = "#26ACB5"), num_features = 15) +
theme_minimal()
valencia.lithics.rf.xv.predictions <- collect_predictions(valencia.lithics.rf.xv.fit, summarize = TRUE) %>%
arrange(.row) %>%
rename(predicted.age = .pred_class,
MP=.pred_MP,
UP=.pred_UP,
EPI=.pred_EPI,
MESO=.pred_MESO,
ENEOL=.pred_ENEOL,
LNEOL=.pred_LNEOL,
true.age=period) %>%
select(-.row, -.config) %>%
relocate(predicted.age, .after = true.age) %>%
bind_cols(training_data_E_Iberia[1], .)
print(valencia.lithics.rf.xv.predictions)
valencia.lithics.rf.xv.predictions %>%
pivot_longer(cols = 4:ncol(valencia.lithics.rf.xv.predictions), names_to = "period", values_to = "probability") %>%
mutate(period = factor(period, levels = c("MP","UP","EPI","MESO","ENEOL","MNEOL","LNEOL"))) %>%
ggplot(aes(x=period, y=probability)) +
geom_line(group=1) +
geom_vline(aes(xintercept = true.age), color="red", size=2, alpha=.5) +
geom_vline(aes(xintercept = predicted.age), color="blue", size=0.8) +
labs(title="Random Forest with Cross-Validated Predictions",
subtitle="Time Periods for Each Assemblage",
x="predicted time period\n(blue line indicates prediction & red line indicates radiocarbon age)") +
facet_wrap(vars(ID)) +
theme_bw(base_size = 16) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
caret::confusionMatrix(
as.factor(valencia.lithics.rf.xv.predictions$true.age),
as.factor(valencia.lithics.rf.xv.predictions$predicted.age))
Confusion Matrix and Statistics
Reference
Prediction MP UP EPI MESO ENEOL LNEOL
MP 51 0 0 0 0 0
UP 0 9 1 0 5 0
EPI 1 3 2 0 0 0
MESO 0 3 0 11 5 0
ENEOL 0 1 0 6 22 0
LNEOL 0 0 0 0 3 6
Overall Statistics
Accuracy : 0.7829
95% CI : (0.7018, 0.8507)
No Information Rate : 0.4031
P-Value [Acc > NIR] : < 2.2e-16
Kappa : 0.7073
Mcnemar's Test P-Value : NA
Statistics by Class:
Class: MP Class: UP Class: EPI Class: MESO Class: ENEOL Class: LNEOL
Sensitivity 0.9808 0.56250 0.66667 0.64706 0.6286 1.00000
Specificity 1.0000 0.94690 0.96825 0.92857 0.9255 0.97561
Pos Pred Value 1.0000 0.60000 0.33333 0.57895 0.7586 0.66667
Neg Pred Value 0.9872 0.93860 0.99187 0.94545 0.8700 1.00000
Prevalence 0.4031 0.12403 0.02326 0.13178 0.2713 0.04651
Detection Rate 0.3953 0.06977 0.01550 0.08527 0.1705 0.04651
Detection Prevalence 0.3953 0.11628 0.04651 0.14729 0.2248 0.06977
Balanced Accuracy 0.9904 0.75470 0.81746 0.78782 0.7771 0.98780
medland.survey.rf.mod <-
rand_forest(trees=500) %>%
set_engine("ranger") %>%
set_mode("classification")
print(medland.survey.rf.mod)
Random Forest Model Specification (classification)
Main Arguments:
trees = 500
Computational engine: ranger
medland.survey.rf.fit <-
medland.survey.rf.mod %>%
fit(as.factor(period) ~ ., data = training_data_E_Iberia[,2:ncol(training_data_E_Iberia)])
print(medland.survey.rf.fit)
parsnip model object
Ranger result
Call:
ranger::ranger(x = maybe_data_frame(x), y = y, num.trees = ~500, num.threads = 1, verbose = FALSE, seed = sample.int(10^5, 1), probability = TRUE)
Type: Probability estimation
Number of trees: 500
Sample size: 129
Number of independent variables: 15
Mtry: 3
Target node size: 10
Variable importance mode: none
Splitrule: gini
OOB prediction error (Brier s.): 0.2341059
medland.survey.rf.predicted <-
medland.survey.rf.fit %>%
predict(medland_survey2014_2017_lithics) %>%
bind_cols(medland_survey2014_2017_info, .,
predict(medland.survey.rf.fit, medland_survey2014_2017_lithics, type="prob")) %>%
rename(predicted.period = .pred_class,
MP=.pred_MP,
UP=.pred_UP,
EPI=.pred_EPI,
MESO=.pred_MESO,
ENEOL=.pred_ENEOL,
LNEOL=.pred_LNEOL) %>%
mutate(predicted.period = factor(predicted.period,
levels=c("MP", "UP", "EPI", "MESO", "ENEOL", "LNEOL")))
print(medland.survey.rf.predicted)
Use the 10th percentile of probabilities from random forest model with dated collections
valencia.lithics.rf.probquantiles <- with(valencia.lithics.rf.xv.predictions %>%
pivot_longer(cols = 4:ncol(valencia.lithics.rf.xv.predictions),
values_to = "probability"),
quantile(probability, probs = c(.1, .25, .5, .75)))
print(paste("10th percentile of random forest probabilities =",
valencia.lithics.rf.probquantiles[1]))
[1] "10th percentile of random forest probabilities = 0.00173557343664358"
# show this in a graph
valencia.lithics.rf.xv.predictions %>%
pivot_longer(cols = 4:ncol(valencia.lithics.rf.xv.predictions), names_to = "period", values_to = "probability") %>%
mutate(period = factor(period, levels = c("MP","UP","EPI","MESO","ENEOL","MNEOL","LNEOL"))) %>%
ggplot(aes(x=probability)) +
geom_density() +
geom_vline(xintercept = valencia.lithics.rf.probquantiles[1], color = 'red') +
geom_vline(xintercept = valencia.lithics.rf.probquantiles[2], color = 'green') +
geom_vline(xintercept = valencia.lithics.rf.probquantiles[3], color = 'blue') +
geom_vline(xintercept = valencia.lithics.rf.probquantiles[4], color = 'green') +
labs(title="Distribution of Combined Random Forest Predictions",
x="age prediction probabilities \n(blue = median, red = 10th percentile, green = 25th and 75th percentile")
# Create base file for graphing and output
medland.survey.rf.graph <-
left_join(select(medland_survey2014_2017, ID, study.area, zone, sector, subsector, total.lithics, area.sqm) %>%
mutate(assemblage = paste(study.area, "-", zone, "-", sector, "-", subsector, sep = "")),
medland.survey.rf.predicted) %>%
dplyr::filter(total.lithics>0) %>%
mutate(density.km2 = total.lithics/area.sqm/1000)
Joining, by = c("ID", "study.area", "assemblage")
# Create output file
# Use 10th percentile for overall ubiquity for patches with only undiagnostic lithics
medland.survey.rf.out <- medland.survey.rf.graph %>%
mutate(MP = replace_na(MP, valencia.lithics.rf.probquantiles[1]),
UP = replace_na(UP, valencia.lithics.rf.probquantiles[1]),
EPI = replace_na(EPI, valencia.lithics.rf.probquantiles[1]),
MESO = replace_na(MESO, valencia.lithics.rf.probquantiles[1]),
ENEOL = replace_na(ENEOL, valencia.lithics.rf.probquantiles[1]),
LNEOL = replace_na(LNEOL, valencia.lithics.rf.probquantiles[1])) %>%
# occupational ubiquity
rename(MP_ubiq = MP,
UP_ubiq = UP,
EPI_ubiq = EPI,
MESO_ubiq = MESO,
ENEOL_ubiq = ENEOL,
LNEOL_ubiq = LNEOL) %>%
# calculate occupational intensity
mutate(MP_int = MP_ubiq*density.km2/800,
UP_int = UP_ubiq*density.km2/250,
EPI_int = EPI_ubiq*density.km2/40,
MESO_int = MESO_ubiq*density.km2/35,
ENEOL_int = ENEOL_ubiq*density.km2/25,
LNEOL_int = LNEOL_ubiq*density.km2/12)
# create file to graph results of patches with diagnostic lithics
medland.survey.rf.graph <- medland.survey.rf.graph %>%
dplyr::filter(!is.na(predicted.period)) %>%
pivot_longer(cols = 10:15,
names_to = "period",
values_to = "ubiquity") %>%
mutate(period = factor(period,
levels = c("MP","UP","EPI","MESO","ENEOL","MNEOL","LNEOL")),
age = case_when(period == "MP" ~ 80000,
period == "UP" ~ 25000,
period == "EPI" ~ 13000,
period == "MESO" ~ 9000,
period == "ENEOL" ~ 6000,
period == "LNEOL" ~ 4000),
predicted.age = case_when(
predicted.period == "MP" ~ 80000,
predicted.period == "UP" ~ 25000,
predicted.period == "EPI" ~ 13000,
predicted.period == "MESO" ~ 9000,
predicted.period == "ENEOL" ~ 6000,
predicted.period == "LNEOL" ~ 4000),
duration.centuries = case_when(
period == "MP" ~ 800,
period == "UP" ~ 250,
period == "EPI" ~ 40,
period == "MESO" ~ 35,
period == "ENEOL" ~ 25,
period == "LNEOL" ~ 12))
medland.survey.rf.graph %>%
ggplot(aes(x=period, y=ubiquity)) +
geom_rect(aes(fill = study.area),xmin = -Inf,xmax = Inf, ymin = -Inf,ymax = Inf,alpha = 0.1) +
geom_line(group=1) +
#geom_vline(aes(xintercept = predicted.period), color="blue") +
labs(title="Random Forest Age Estmates for Medland Survey Data",
subtitle="Time Periods for Each Assemblage, Colored by Study Area",
x="predicted time period\n(blue line indicates maximum predicted probability)") +
facet_wrap(~ assemblage + ID) +
theme_bw(base_size = 16) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
medland.survey.rf.graph %>%
ggplot(aes(x=period, y=ubiquity*density.km2/duration.centuries)) +
geom_rect(aes(fill = study.area),xmin = -Inf,xmax = Inf, ymin = -Inf,ymax = Inf,alpha = 0.1) +
geom_line(group=1) +
#geom_vline(aes(xintercept = predicted.period), color="blue") +
labs(title="Land Use Intensity Estmates for Medland Survey Data",
subtitle="Time Periods for Each Assemblage, Colored by Study Area",
x="predicted time period\n(blue line indicates maximum predicted probability)",
y="estimated artifact accumulation rate\n(artifacts/km2/century)") +
facet_wrap(~ assemblage + ID) +
theme_bw(base_size = 16) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
medland.survey.rf.graph %>%
dplyr::filter(ID==238 | ID==960) %>%
ggplot(aes(x=period, y=ubiquity)) +
geom_line(group=1, size=2) +
labs(title="Occupational Ubiquity for 2 Survey Patches",
subtitle="Random Forest Age Estmates",
x="predicted time period",
y="occupational ubiquity") +
facet_wrap(~ assemblage) +
theme_bw(base_size = 20) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
medland.survey.rf.graph %>%
dplyr::filter(ID==238 | ID==960) %>%
ggplot(aes(x=period, y=ubiquity*density.km2/duration.centuries)) +
geom_line(group=1, size=2) +
labs(title="Land Use Intensity for 2 Survey Patches",
subtitle="Random Forest Age Estmates",
x="predicted time period",
y="artifacts / km^2 /century") +
facet_wrap(~ assemblage) +
theme_bw(base_size = 20) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
axis.title.y = element_markdown())
write_csv(medland.survey.rf.out, "medland_survey_rf.csv")
totals.surveyed <- medland_survey2014_2017 %>%
select(study.area, area.sqm) %>%
group_by(study.area) %>%
summarise(total.area.sqm=sum(area.sqm))
left_join(medland.survey.rf.graph, totals.surveyed) %>%
dplyr::filter(ubiquity > quantile(medland.survey.rf.graph$ubiquity,
probs = .5)) %>%
group_by(study.area, period) %>%
summarize(area.sum = sum(area.sqm),
total.area.sqm = first(total.area.sqm)) %>%
select(study.area, period, area.sum, total.area.sqm) %>%
mutate(pct.coverage = area.sum/total.area.sqm) %>%
ungroup() %>%
add_row(study.area = "Canal de Navarrés",
period = "EPI",
area.sum = 0,
total.area.sqm = 4107582,
pct.coverage = 0) %>%
add_row(study.area = "Cocina/Catadau",
period = "EPI",
area.sum = 0,
total.area.sqm = 1095683,
pct.coverage = 0) %>%
add_row(study.area = "Cocina/Catadau",
period = "MESO",
area.sum = 0,
total.area.sqm = 1095683,
pct.coverage = 0) %>%
add_row(study.area = "Hoya de Buñol",
period = "EPI",
area.sum = 0,
total.area.sqm = 1551377,
pct.coverage = 0) %>%
# Needed for upper quartile graphing
# add_row(study.area = "Canal de Navarrés",
# period = "MP",
# area.sum = 0,
# total.area.sqm = 4107582,
# pct.coverage = 0) %>%
# add_row(study.area = "Canal de Navarrés",
# period = "MESO",
# area.sum = 0,
# total.area.sqm = 4107582,
# pct.coverage = 0) %>%
# add_row(study.area = "Canal de Navarrés",
# period = "LNEOL",
# area.sum = 0,
# total.area.sqm = 4107582,
# pct.coverage = 0) %>%
# add_row(study.area = "Cocina/Catadau",
# period = "MP",
# area.sum = 0,
# total.area.sqm = 1095683,
# pct.coverage = 0) %>%
# add_row(study.area = "Cocina/Catadau",
# period = "LNEOL",
# area.sum = 0,
# total.area.sqm = 1095683,
# pct.coverage = 0) %>%
# add_row(study.area = "Hoya de Buñol",
# period = "MP",
# area.sum = 0,
# total.area.sqm = 1551377,
# pct.coverage = 0) %>%
# add_row(study.area = "Hoya de Buñol",
# period = "MESO",
# area.sum = 0,
# total.area.sqm = 1551377,
# pct.coverage = 0) %>%
# add_row(study.area = "Hoya de Buñol",
# period = "LNEOL",
# area.sum = 0,
# total.area.sqm = 1551377,
# pct.coverage = 0) %>%
ggplot(aes(x=factor(period, levels =
c('MP','UP','EPI','MESO','ENEOL','LNEOL')),
y=pct.coverage, group=1)) +
geom_line(lwd=1.5) +
facet_wrap(~study.area, ncol = 1) +
labs(title='Aggregate Occupational Ubiquity',
subtitle='Proportion of study area surveyed\nwith ubiquity above the median',
x='period',
y='proportion') +
theme_bw(base_size = 16)
Joining, by = "study.area"`summarise()` has grouped output by 'study.area'. You can override using the `.groups` argument.
left_join(medland.survey.rf.graph, totals.surveyed) %>%
mutate(intensity = ubiquity*density.km2/duration.centuries) %>%
dplyr::filter(intensity >= quantile(medland.survey.rf.graph$ubiquity*medland.survey.rf.graph$density.km2/medland.survey.rf.graph$duration.centuries, probs =.5)) %>%
group_by(study.area, period) %>%
summarize(area.sum = sum(area.sqm),
total.area.sqm = first(total.area.sqm),
density.km2 = first(density.km2)) %>%
select(study.area, period, area.sum, total.area.sqm) %>%
mutate(pct.coverage = area.sum/total.area.sqm) %>%
ungroup %>%
add_row(study.area = "Canal de Navarrés",
period = "MP",
area.sum = 0,
total.area.sqm = 4107582,
pct.coverage = 0) %>%
add_row(study.area = "Hoya de Buñol",
period = "MP",
area.sum = 0,
total.area.sqm = 1551377,
pct.coverage = 0) %>%
# Needed for optional upper quartile graphing
# add_row(study.area = "Hoya de Buñol",
# period = "UP",
# area.sum = 0,
# total.area.sqm = 1551377,
# pct.coverage = 0) %>%
ggplot(aes(x=factor(period, levels =
c('MP','UP','EPI','MESO','ENEOL','LNEOL')),
y=pct.coverage, group=1)) +
geom_line(lwd=1.5) +
facet_wrap(~study.area, ncol = 1) +
labs(title='Aggregate Land Use Intensity',
subtitle='Proportion of study area surveyed\nwith intensity above the median',
x='period',
y='proportion') +
theme_bw(base_size = 16)
Joining, by = "study.area"`summarise()` has grouped output by 'study.area'. You can override using the `.groups` argument.
Only use dates with COV ≤ 0.05
C14_SE_Iberia_all <- C14_SE_Iberia_all %>% dplyr::filter(C14.CV<0.05 & C14.SD>0)
all.dates.calibrated <- with(C14_SE_Iberia_all, BchronCalibrate(ages = C14.mean, ageSds = C14.SD, calCurves = calib.curve, positions = site))
C14_SE_Iberia_all$BP.cal.median <- sapply(1:length(all.dates.calibrated), function(x) round(median(all.dates.calibrated[[x]]$ageGrid)))
C14_SE_Iberia_all.bins <- C14_SE_Iberia_all %>%
with(., binPrep(site, C14.mean, 100))
C14_SE_Iberia_all.modeltest <- C14_SE_Iberia_all %>%
with(., calibrate(x=C14.mean, errors=C14.SD, calCurves = calib.curve, normalised=TRUE, calMatrix=FALSE)) %>%
modelTest(.,
errors = C14_SE_Iberia_all$C14.SD,
timeRange = c(35000,3000),
runm = 500,
model="exponential",
datenormalised=TRUE,
nsim = 200,
ncores = ncores,
method = 'calsample',
bins = C14_SE_Iberia_all.bins)
[1] "Calibrating radiocarbon ages..."
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 3%
|
|=== | 4%
|
|==== | 4%
|
|==== | 5%
|
|===== | 5%
|
|===== | 6%
|
|====== | 6%
|
|====== | 7%
|
|======= | 7%
|
|======= | 8%
|
|======== | 8%
|
|======== | 9%
|
|========= | 9%
|
|========= | 10%
|
|========== | 10%
|
|========== | 11%
|
|=========== | 11%
|
|=========== | 12%
|
|============ | 12%
|
|============ | 13%
|
|============= | 13%
|
|============= | 14%
|
|============== | 15%
|
|============== | 16%
|
|=============== | 16%
|
|=============== | 17%
|
|================ | 17%
|
|================ | 18%
|
|================= | 18%
|
|================= | 19%
|
|================== | 19%
|
|================== | 20%
|
|=================== | 20%
|
|=================== | 21%
|
|==================== | 21%
|
|==================== | 22%
|
|===================== | 22%
|
|===================== | 23%
|
|====================== | 23%
|
|====================== | 24%
|
|======================= | 24%
|
|======================= | 25%
|
|======================== | 25%
|
|======================== | 26%
|
|========================= | 26%
|
|========================= | 27%
|
|========================== | 27%
|
|========================== | 28%
|
|=========================== | 29%
|
|=========================== | 30%
|
|============================ | 30%
|
|============================ | 31%
|
|============================= | 31%
|
|============================= | 32%
|
|============================== | 32%
|
|============================== | 33%
|
|=============================== | 33%
|
|=============================== | 34%
|
|================================ | 34%
|
|================================ | 35%
|
|================================= | 35%
|
|================================= | 36%
|
|================================== | 36%
|
|================================== | 37%
|
|=================================== | 37%
|
|=================================== | 38%
|
|==================================== | 38%
|
|==================================== | 39%
|
|===================================== | 39%
|
|===================================== | 40%
|
|====================================== | 40%
|
|====================================== | 41%
|
|======================================= | 41%
|
|======================================= | 42%
|
|======================================== | 43%
|
|======================================== | 44%
|
|========================================= | 44%
|
|========================================= | 45%
|
|========================================== | 45%
|
|========================================== | 46%
|
|=========================================== | 46%
|
|=========================================== | 47%
|
|============================================ | 47%
|
|============================================ | 48%
|
|============================================= | 48%
|
|============================================= | 49%
|
|============================================== | 49%
|
|============================================== | 50%
|
|=============================================== | 50%
|
|=============================================== | 51%
|
|================================================ | 51%
|
|================================================ | 52%
|
|================================================= | 52%
|
|================================================= | 53%
|
|================================================== | 53%
|
|================================================== | 54%
|
|=================================================== | 54%
|
|=================================================== | 55%
|
|==================================================== | 55%
|
|==================================================== | 56%
|
|===================================================== | 56%
|
|===================================================== | 57%
|
|====================================================== | 58%
|
|====================================================== | 59%
|
|======================================================= | 59%
|
|======================================================= | 60%
|
|======================================================== | 60%
|
|======================================================== | 61%
|
|========================================================= | 61%
|
|========================================================= | 62%
|
|========================================================== | 62%
|
|========================================================== | 63%
|
|=========================================================== | 63%
|
|=========================================================== | 64%
|
|============================================================ | 64%
|
|============================================================ | 65%
|
|============================================================= | 65%
|
|============================================================= | 66%
|
|============================================================== | 66%
|
|============================================================== | 67%
|
|=============================================================== | 67%
|
|=============================================================== | 68%
|
|================================================================ | 68%
|
|================================================================ | 69%
|
|================================================================= | 69%
|
|================================================================= | 70%
|
|================================================================== | 70%
|
|================================================================== | 71%
|
|=================================================================== | 72%
|
|=================================================================== | 73%
|
|==================================================================== | 73%
|
|==================================================================== | 74%
|
|===================================================================== | 74%
|
|===================================================================== | 75%
|
|====================================================================== | 75%
|
|====================================================================== | 76%
|
|======================================================================= | 76%
|
|======================================================================= | 77%
|
|======================================================================== | 77%
|
|======================================================================== | 78%
|
|========================================================================= | 78%
|
|========================================================================= | 79%
|
|========================================================================== | 79%
|
|========================================================================== | 80%
|
|=========================================================================== | 80%
|
|=========================================================================== | 81%
|
|============================================================================ | 81%
|
|============================================================================ | 82%
|
|============================================================================= | 82%
|
|============================================================================= | 83%
|
|============================================================================== | 83%
|
|============================================================================== | 84%
|
|=============================================================================== | 84%
|
|=============================================================================== | 85%
|
|================================================================================ | 86%
|
|================================================================================ | 87%
|
|================================================================================= | 87%
|
|================================================================================= | 88%
|
|================================================================================== | 88%
|
|================================================================================== | 89%
|
|=================================================================================== | 89%
|
|=================================================================================== | 90%
|
|==================================================================================== | 90%
|
|==================================================================================== | 91%
|
|===================================================================================== | 91%
|
|===================================================================================== | 92%
|
|====================================================================================== | 92%
|
|====================================================================================== | 93%
|
|======================================================================================= | 93%
|
|======================================================================================= | 94%
|
|======================================================================================== | 94%
|
|======================================================================================== | 95%
|
|========================================================================================= | 95%
|
|========================================================================================= | 96%
|
|========================================================================================== | 96%
|
|========================================================================================== | 97%
|
|=========================================================================================== | 97%
|
|=========================================================================================== | 98%
|
|============================================================================================ | 98%
|
|============================================================================================ | 99%
|
|=============================================================================================| 100%
[1] "Done."
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6 ✔ purrr 0.3.4
✔ tibble 3.1.7 ✔ dplyr 1.0.9
✔ tidyr 1.2.0 ✔ stringr 1.4.0
✔ readr 2.1.2 ✔ forcats 0.5.1
Attaching package: ‘snow’
The following objects are masked from ‘package:parallel’:
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
parCapply, parLapply, parRapply, parSapply, splitIndices,
stopCluster
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6 ✔ purrr 0.3.4
✔ tibble 3.1.7 ✔ dplyr 1.0.9
✔ tidyr 1.2.0 ✔ stringr 1.4.0
✔ readr 2.1.2 ✔ forcats 0.5.1
Attaching package: ‘snow’
The following objects are masked from ‘package:parallel’:
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
parCapply, parLapply, parRapply, parSapply, splitIndices,
stopCluster
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6 ✔ purrr 0.3.4
✔ tibble 3.1.7 ✔ dplyr 1.0.9
✔ tidyr 1.2.0 ✔ stringr 1.4.0
✔ readr 2.1.2 ✔ forcats 0.5.1
Attaching package: ‘snow’
The following objects are masked from ‘package:parallel’:
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
parCapply, parLapply, parRapply, parSapply, splitIndices,
stopCluster
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6 ✔ purrr 0.3.4
✔ tibble 3.1.7 ✔ dplyr 1.0.9
✔ tidyr 1.2.0 ✔ stringr 1.4.0
✔ readr 2.1.2 ✔ forcats 0.5.1
Attaching package: ‘snow’
The following objects are masked from ‘package:parallel’:
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
parCapply, parLapply, parRapply, parSapply, splitIndices,
stopCluster
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6 ✔ purrr 0.3.4
✔ tibble 3.1.7 ✔ dplyr 1.0.9
✔ tidyr 1.2.0 ✔ stringr 1.4.0
✔ readr 2.1.2 ✔ forcats 0.5.1
Attaching package: ‘snow’
The following objects are masked from ‘package:parallel’:
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
parCapply, parLapply, parRapply, parSapply, splitIndices,
stopCluster
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6 ✔ purrr 0.3.4
✔ tibble 3.1.7 ✔ dplyr 1.0.9
✔ tidyr 1.2.0 ✔ stringr 1.4.0
✔ readr 2.1.2 ✔ forcats 0.5.1
Attaching package: ‘snow’
The following objects are masked from ‘package:parallel’:
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
parCapply, parLapply, parRapply, parSapply, splitIndices,
stopCluster
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6 ✔ purrr 0.3.4
✔ tibble 3.1.7 ✔ dplyr 1.0.9
✔ tidyr 1.2.0 ✔ stringr 1.4.0
✔ readr 2.1.2 ✔ forcats 0.5.1
Attaching package: ‘snow’
The following objects are masked from ‘package:parallel’:
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
parCapply, parLapply, parRapply, parSapply, splitIndices,
stopCluster
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
✔ ggplot2 3.3.6 ✔ purrr 0.3.4
✔ tibble 3.1.7 ✔ dplyr 1.0.9
✔ tidyr 1.2.0 ✔ stringr 1.4.0
✔ readr 2.1.2 ✔ forcats 0.5.1
Attaching package: ‘snow’
The following objects are masked from ‘package:parallel’:
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, clusterSplit, makeCluster, parApply,
parCapply, parLapply, parRapply, parSapply, splitIndices,
stopCluster
[1] "Aggregating observed dates..."
[1] "Monte-Carlo test..."
par(mar=c(7,7,7,3))
plot(C14_SE_Iberia_all.modeltest, xlim = c(30000,4000), col.obs = 'black', lwd.obs = 5, drawaxes = F)
axis(1, cex.axis = 2, pos = -.01, at=(seq(30000,0, by=-5000)))
axis(2, cex.axis = 2, pos = 30200)
mtext(side=1, line=5, "calibrated years BP", cex=2.5)
mtext(side=2, line=4, "summed probability density", cex=2.5)
title(main=paste("Southern and Eastern Iberia SPD (N = ", nrow(C14_SE_Iberia_all), ")\n"), cex.main = 3)